Two-dimensional flagellar synchronization in viscoelastic fluids 
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Experimental studies have demonstrated that spermatozoa synchronize their flagella when swim- 
ming in close proximity. In a Newtonian fluid, it was shown theoretically that such synchronization 
arises passively due to hydrodynamic forces between the two swimmers if their waveforms exhibit 
a front-back geometrical asymmetry. Motivated by the fact that most biological fluids possess a 
polymeric microstructure, we address here synchronization in a viscoelastic fluid analytically. Using 
a two-dimensional infinite sheet model we show that the presence of polymeric stresses removes 
the geometrical asymmetry constraint, and therefore even symmetric swimmers synchronize. Such 
synchronization occurs on asymptotically faster time scales than in a Newtonian fluid, and the swim- 
mers are seen to be driven into a stable in-phase conformation minimizing the energy dissipated in 
the surrounding fluid. 

INTRODUCTION 

Swimming microorganisms are found everywhere in nature, from bacteria moving towards nutrients, to spermatozoa 
migrating towards an ovum. Given their size, they swim in a realm where viscous stresses predominate and inertia is 
negligible. This fact often necessitates strategies for locomotion which are vastly different from those of organisms in 
the macroscopic world [1]. Much progress has been made in the physical and hydrodynamic analysis of microscopic 
swimming [2-4], yet many marked phenomena are still of current interest. One such phenomenon is the experimentally- 
observed synchronization of spermatozoa flagella (slender filaments deformed in a wavelike fashion and propelling the 
cells forward) when swimming in close proximity [5, 6]. This synchronization has been observed to prompt an 
increase in the swimming speed of the co-swimming cells which could thereby provide a competitive, and hence 
perhaps evolutionary, advantage [6] . The theoretical analysis of the synchronization of swimming microorganisms in 
Newtonian fluids dates back to the work of Taylor who, assuming two infinite and parallel two-dimensional sheets, 
showed that synchronous beating minimizes the viscous dissipation between the model cells [7] . This synchronization 
was subsequently demonstrated computationally for both infinite [8] , and finite two-dimensional models [9] . Recently, 
it was shown that phase-locking arises purely passively due to the fluid forces between the two swimmers, and requires 
a front-back asymmetry in the flagellar geometry of the cells without which no synchronization can occur [10]. 

In all previously-studied situations, synchronization was addressed in the case of a Newtonian fluid. However, 
most biological fluids involved, for example, in mammalian reproduction are non-Newtonian. As mammalian sper- 
matozoa make their journey through the female reproductive tract they encounter several complex fluids, including 
glycoprotein-based cervical mucus in the cervix, mucosal epithelium inside the fallopian tubes, and actin-based vis- 
coelastic gel outside the ovum [11, 12]. In a viscoelastic fluid, kinematic reversibility, restated in Purcell's scallop 
theorem [1], breaks down due to the presence of normal stresses and shear-dependent material functions, fundamen- 
tally altering the governing flow physics [13]. The waveform, structure, and swimming path of spermatozoa have been 
experimentally observed to be modified in viscoelastic fluids [14]. Locomotion in complex fluids has been studied 
analytically [15-17], and it has been shown that microorganisms which propel themselves by propagating waves along 
their flagella have a lower swimming speed in a viscoelastic fluid than in a Newtonian fluid [18, 19]. 

In this paper we study the passive synchronization of two flagellated cells in a viscoelastic (Oldroyd-B) fluid. Using 
Taylor's infinite two-dimensional sheet model, we show that not only does phase locking arise in a viscoelastic fluid, 
but also that it does not require the front-back geometrical asymmetry that must exist for such a model to display 
synchronization in a Newtonian fluid. We demonstrate that the system evolves to a single stable fixed point at the 
in-phase conformation, which is also the conformation that yields minimal energy dissipation. In addition, we show 
that the evolution to a phase-locked state occurs on asymptotically faster time scales than in a Newtonian fluid. 
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FIG. 1: Model system consisting of two infinite sinusoidal sheets passing waves at speed c and thereby swimming at speed 
U±Ua/2 in the opposite direction. The difference in phase <f> incurs a relative velocity Ua(0, h) between the two sheets denoted 
positive when the top sheet swims to the right relative to the bottom one. The average separation distance is denoted h. 



SETUP 



Our system, shown in Fig. 1, consists of two parallel infinite two-dimensional sheets, separated by a mean distance 
h. Both sheets propagate sinusoidal waves of transverse displacement of amplitude a at speed c = uj/k, where w 
is the wave frequency and k is the wavenumber, but have an initial phase difference <f>Q. By passing these waves, 
the sheets propel themselves in the direction opposite to the wave speed [7]. The sheets are also permitted to move 
relative to each other with an unknown velocity U A , denoted positive when the top sheet (#2) swims in the positive 
x direction relative to the bottom one (#1). The positions of the sheets, in their swimming frames, are thereby given 
by y x = asm(kx - uit - (j> /2 + J* kU A (t')dt' /2) and y 2 = h + asm(kx -ojt+ (j> /2 - J* kU A (t')dt' /2). 

We use the following dimcnsionless variables x* — xk, t* — tui, u* = u/c, U A = U A /c. The amplitude of the waves 
is non-dimensionalized by the wavenumber, e = ak. For convenience we let x* = x* — t* and <f> = (f>o — f Q U A (t')dt' 
which is the instantaneous phase difference between the two sheets. The position of the sheets in dimensionless form is 
thus given by y\ = e sin(x* — <fi/2), and y| = h* + e sin(x* + cf>/2), and the phase evolves in time according to <f> = —U A . 
We refer to the <f> = conformation as in-phase, and the <p = n conformation as opposite-phase. The system is 2-7T 
periodic and (f> is defined from — tt to tt. We now drop the (*) notation and refer only below to dimcnsionless variables. 

Since this problem is two dimensional, we introduce the streamfunction ip(x,y,t) where the components of the 
velocity field are u = [u, v] T — [dt/j/dy,~dip/dx] T , and the incompressibility condition is always satisfied. The 
boundary conditions become then 



V4>\ y=yi = V -U A y/2 + esin{x-(f>/2) , (1) 

L J y=vi 

ViP\ y=y2 = v\u A y/2 + esm{x + (j)/2)\ . (2) 

Mechanical equilibrium for a low-Reynolds number flow in a viscoelastic fluid is given by Vp = V • t, where p is 
the pressure and t is the deviatoric part of the stress tensor. To relate the stress to the strain-rate, we use in this 
paper the simplest polymeric constitutive equation derived from a microscopic theory, namely the Oldroyd-B model 
[20, 21] given by 



v 



V 



r + Dei t= 7 + De 2 7, (3) 

where 7 = Vu + Vu T is the strain-rate tensor; the upper convected derivative is defined for a general tensor A as 
A= dA/dt + u • VA — (Vu T • A + A • Vu). We have defined two Deborah numbers, Dei = Aiu;, where Ai is the 
relaxation time of the polymer, and De 2 = where A2 = Xir) s /(7] s + rj p ) is the retardation time of the polymer and 
•q Sl rj p refer to the contribution to the viscosity from the solvent and polymer respectively [21]. Rheological studies 
have shown the relaxation time of cervical mucus to be between 1 and 10 seconds [22]. Given that the flagella of 
spermatozoa typically beat at a frequency between 20 and 50 Hz [3] we might expect a range of Deborah numbers 
Dei = 10 2 — 10 3 , and in most practical instances Dei 3> De 2 [21]. 

In the spirit of Taylor's seminal study, we look to solve this system in a small amplitude limit, e <C 1, with a 
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perturbation expansion in powers of e for stress r = et 1 + e 2 T2 + strain-rate 7 = ej x + e 2 j 2 + streamfunction 
ip = etpi + e 2 ij)2 + ... and the relative velocity between the sheets U\ = cUai + £ 2 Ua2 + ■■■■ 

ANALYSIS 
First-order solution 



The leading order component of (3) is 



dTi 

dt 



Ti + Dci^^- = j 1 + Dc 2 



dt ' 



Taking the divergence and the curl of (4) we get the governing equation for the first-order streamfunction 



1 + De 2 ^ ) V J n = 0. 



With the first order boundary conditions 



VV>i| w =o = V[(-£/'Ai»/2 + sin(a:-^/2)]| w= o, 
^i\y=h = V[U A iy/2 + sm{x + (t>/2)]\ y = h , 



(4) 



(5) 



(6) 
(7) 



the solution can be shown to be 



ipi = ao{y) + ai(y) cos(x) + b x {y) sin(x), 



where 



ao(y) = C iy 2 (y - h\ + l -U Al y (| - l) , 



ai(y) = 



h(y) 



h — sinh(/i) 



cos ■ 



h + sinh(/i) 



/ cosh(/i — y) — (h — y) cosh y + sinh(/i — y) — sinh y^j , 
(y cosh(/i — y) + (h — y) cosh y + sinh(/i — y) + sinh y^j . 



(8) 

(9) 
(10) 
(11) 



To determine the unknown constant Ci we resort to dynamical considerations. For simplicity we resolve the 



^ 0) + ^ 



where 



streamfunction into its complex Fourier components in the wave variable x = x — t giving ipi = 5R 

denotes the real part and ipi^ — a (y) and Tp{^ = {ai{y) J rib\{y))e~ lx . The strain-rate tensor, A f 1 = Vui+Vu^, 
can then be obtained using (8). Exploiting (4), we see that the first-order stress tensor is given by 



n = 5ft 



, (o) l + tDea wi) 
71 + l + zD Cl 71 



(12) 



If <j = —pi + t refers to the total stress tensor, integration of V • er = leads to the sum of the forces, f , on the upper 
and lower sheets (over a period) equal to zero, i.e. f\ y = yi + f\ y = y2 — 0. At leading order, the horizontal component 
of this relationship is 



/•27T />Z7T 
/ Tl X y\y = dx = / Tl X y\y=hdX, (13) 

Jo Jo 

which yields Ci = 0. We finally determine the relative velocity by insisting the sheets be force-free. Typically for 
each sheet, one must sum the forces on both the inner and outer surfaces. However, the outer problem is force-free 
for all Uai [18]. The net force on the upper sheet is therefore given by 



r 2n 



r2n 

fix = - / a' \h)dx = -2nU A1 /h, 
Jo 



(14) 
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and hence Uai is zero, which is expected due to the e — > — e symmetry of the system. With Ci and J7ai equal to 
zero then ao — (see (9)) therefore we have no time-averaged flow and we get a simplified relation between stress and 
strain-rate in Fourier space as 



f + iBe 2 r 

Tl TTTdcT 71 ' 



(15) 



Second-order solution 



The second-order component of (3) is given by 

(\ + Dei T 2 -(l + De 2 7 2 = Dei (V< • n + n ■ Vui - m • Vr) 

- Dc 2 (Vu^ ■ 7! + 7! ■ Vui - ui • V-Vi) . 



(16) 



The only part of the streamfunction, ^2, that will contribute to the force on the sheets at second order is its mean 
value in x. Using (15), the mean value of (16) is given by 



<T2) - <7 2 > = » 



Dei-De 2 ,__ T 



2(l + iDei) 



(VuJ* • 7i + 7i ' VuJ - u{ ■ V 7l ) 



(17) 



where *'s indicate complex conjugates, and (...) denotes averaging over one period in x. The right hand side of (17) 
can then be computed using the first-order streamfunction. Upon taking the divergence and the curl of (17), we 
obtain 



V74/ , \ D °i - De 2 d2 m u 

v<^H TT ^ r ^G(^, 



(18) 



where 



G(y;h,t 



- ai(y)(Deiai'(i/) + 3b'{(y)) + &'i0/) ^'/(j/) - Dei&i'fo) 
+ a 1 (»)(-26' 1 (l/)-De 1 ai"(») + 6i"(i/)) 
+ 6i(l/)(2a' 1 (»)-ai"(i/)-Dei6i"(i,)) 
The second-order components of the boundary conditions are 



Vip 2 \ y =Q = -V(U A2 y)-sm{x-(t)/2)V 



dy 



V^| w=fc = V(U A2 y) - sm(x + 0/2) V (jj^ 



y=h 



Taking the mean value of (20) and (21) yields 

dfo) 



dx 



dy 

9(V>2> 



dx 



y=h 



dy 



y=h 



= o, 
1 

2 

= 0, 

1 

2 



cos 2 (f )(-/i + sinh ft) sin 2 (f)(h + sinh/i) 
-£^A2 H ; — ; — r~r~, H 



/i + sinh ft 



— ft + sinh ft 



t^A2 + 



cos 2 ( f ) (- ft + sinh ft) sin 2 ( § ) (/i + sinh ft) 



ft + sinh ft 



+ 



-ft + sinh ft 



(19) 

(20) 
(21) 

(22) 
(23) 

(24) 
(25) 
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Solving (18) with the above boundary conditions leads to the solution 
<V>2> = C 2 y 2 (y-~— 



3ft \ U A2 y(y-h) i/cos 2 (f)(-ft + sinhft) y sin 2 (f)(ft + sinh ft) 



2ft 



+ 



2(ft + sinhft) 



+ 



2(-ft + sinh ft) 



+ 



Dei - De 2 
1 + Dc 2 



y(y - 2ft) 

2ft 



I Gdy\ y=0 - |^ y Gdi/| w=fc + || Gdj/ S 



(26) 



To find the unknown constant C 2 we again turn to dynamical considerations. Using integration by parts, it is 
straightforward to get that the force on the bottom sheet, to 0(e 2 ), is given by 



2.T 



(T2xy)\y=odx, 



(27) 



and only the mean component of the second-order stress, (t2 XV ), contributes to the net force. A similar relationship 
holds for the force on the upper sheet. We then proceed by obtaining (T 2x y) from (17), where ( / y 2 ) = V(u 2 ) + V(u 2 ) T 



and (u 2 ) = [d(tlj2)/dy,—d(ip2)/dx] T . Exploiting that f\ y = yi 
upper sheet is finally given by 



-f |v=«2 we obtain C 2 = 0, and the net force on the 



2ttC/a2 



2:t 



■ 47T 



Dei — De 2 
1 + De 2 



A(h) sin <p, A(h) 



ft cosh ft + sinh ft 
cosh(2ft) - 2ft 2 - ] 



(28) 



RESULTS 
Synchronization 

It is insightful to first consider the nature of the force which arises if the sheets are not permitted to move relative to 
each other but instead held with a fixed phase difference. If 1/^2 = 0, then the force in (28) is zero for = 0, ±7r. The 
function A(h), governing the variation in the force amplitude with mean distance ft, is positive definite and decays 
exponentially with ft, while becoming unbounded near ft = (see Fig. 2a). Since Dei > De 2 , we see that the force 
f 2x oc sin0. This indicates that (f> — is a stable fixed point while <j> = ±7r are unstable, and therefore we expect 
in-phase synchronization to occur. 

We next observe that we obtain here a nonzero force on sheets with front/back symmetric waveforms [24]. In 
the case of a Newtonian fluid, this is forbidden because of kinematic reversibility, and the force is identically zero 
unless the front/back symmetry is broken [10]. Indeed, in a Newtonian fluid, for any system with both vertical and 
horizontal symmetry one can reflect about both axes of symmetry then reverse the kinematics to obtain an identical 
conformation with the opposite force necessitating f x = (our calculations confirm this by setting Dei = De 2 = in 
28). In a viscoelastic fluid, time is no longer merely a parameter, and therefore the flow is no longer kinematically 
reversible, thereby permitting a nonzero force. 

If instead of holding the sheets fixed, we let them move, we then have to enforce the force-free condition, and thus 
we obtain the relative speed Ua2 = ft/fx/^ 71 "; where f 2x is the static force incurred when Ua2 = in (28). The 
remarkable result is that, since the force occurs in a viscoelastic fluid at second order in the wave amplitude, the 
phase will evolve on a time scale varying as t ~ e~ 2 . In a Newtonian fluid, it can be shown that the force is always 
zero to second order in e, for any shape, and first appears at fourth order for shapes with broken front-back symmetry 
[10]. This means that in a Newtonian fluid, at best, the phase will evolve to a phase-locked configuration on a time 
scale varying as t ~ e~ 4 . In complex fluids, synchronization is therefore seen to take place on asymptotically faster 
time scales than in a Newtonian fluid. 

We now solve analytically for the time-evolution of the phase. Since to leading order </> = — £ 2 C/a2, we obtain a 
differential equation for the evolution for 4> as 

(29) 

which, for constant ft, can be integrated to yield an analytical formula for the phase as 

^2hA(h) f DCl ,^ D 2 C2 V 



</>(t) = 2 tan 1 ^ tan ( ) exp 



1 + De^ 



(30) 
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FIG. 2: a) Amplitude A of the phase- locking force decays exponentially with the separation distance, h (28). b) Time-evolution 
of the phase angle (j>(t) from various initial conditions towards the stable in-phase conformation ((30) with e = 0.1, h — 2, 
Dei = 100 and De 2 = 10). 



Given that Dei > De 2 , we see that cj> ~ ±e~* near the <f> = fixed point, meaning it is stable; however, near the 
<f> = ±7r fixed points, (/> =p n ~ =pe meaning they are unstable and hence the phase converges to = for all initial 
conditions. The time-evolution of the phase from various initial positions, assuming a constant separation between 
the sheets of h = 2, is plotted in Fig. 2b with Dei = 100, De 2 = 10 and e — 0.1. All initial conformations evolve to 
stable in-phase synchrony. In the Newtonian case, the stability of the in-phase versus opposite-phase conformation is 
purely a matter of geometry, regardless of considerations of energy dissipation, and in fact two swimmers can evolve 
to a stable conformation which maximizes the energy dissipated [10]. In contrast, in a viscoelastic fluid we find that 
with no asymmetry the system naturally evolves to an in-phase conformation which, as we show below, coincides with 
the conformation of minimal viscous dissipation. 



Energy dissipation 



The energy dissipated in the fluid between two sheets is given by integrating the dissipation density, r : 7, over the 
volume. The leading order component is given by 

n : 7, = It [ft] : R m = : 7, - ft ' t J . ( 31 ) 

where 9f[...] denotes the imaginary part. The second term in (31) integrates to zero over a period, thus, to leading 
order, the energy dissipation rate per unit depth over one period is given by 



2tt py 2 



^ ^l + De^ / / -V.^i' 1 " 1 "- <*->> 

The result of (32) is merely a scalar multiple of the Newtonian dissipation calculated by [7], which is minimum at the 
in-phase conformation, and maximum in the case of opposite-phase (and decays to zero as h — ► 0). We see therefore 
that in a viscoelastic fluid the system is driven towards a state of minimum energy dissipation. 



Vertical force 



Since the evolution of the phase depends on the separation distance h, it is informative to analyze the magnitude 
of vertical forces between the sheets. We now proceed to compute the vertical force from the inner problem with 
the first and second-order streamfunctions derived here, and we use the solutions of the outer flow problem from the 
literature [18]. The vertical force on the bottom sheet to first order is given by 

fi y = cri 22 | a= o - o-i 2 u 2 tcl | y=0 da;. (33) 
Jo 
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FIG. 3: a) The leading-order vertical force on the lower sheet, f : 2 y (h), displays a complex behavior which depends on both 
phase difference <j> and mean separation h. Plotted for <f> = {0, 7r/16, 7r/8, 7r/4, tv/2 (shown dashed), 37r/4, 7r} (arbitrary units), 
b) Equilibrium separation, heq, defined as the distance at which the vertical force is zero for a given </>. 



Both components are individually zero, hence the force is zero. At second order, the outer flow yields no force for all 
Ua2 therefore the force on the bottom sheet is given by 



hy 



( T 2yy) 



9(T2yy) 

dy 



dy 



dx = 2tt 



y=o 



Dei — De2 
I + De 2 



Bi(h)+B 2 (h) cos ( 



where 



B 1 (h) 
B 2 (h) 



2 [1 - (1 + 2h 2 ) cosh(2/i) - 2/isinh(2/i)] 

(cosh(2/i) - 1 - 2h 2 ) 2 
(Ah 2 - 1) cosh h + cosh(3/i) + 2h [3 + 2h 2 + cosh(2ft,)] sinh h 
(cosh(2/i) - 2h 2 - l) 2 



(34) 

(35) 
(36) 



The vertical component of the force is a cosine function in <j) which is amplified by the positive-definite function 
B 2 and shifted by the negative-definite function B\. Both functions become unbounded as h — ► 0, and both tend 
asymptotically to zero as h — > 00. In Fig. 3a we plot the vertical force as a function of the distance between the 
sheets, h, for various <fi (arbitrary units). If the phase difference is above 7r/2 (<fr — tt/2 is shown dashed), then the 
sheets will be repelled from each other. However, as the sheets get closer in phase there arises a finite equilibrium 
separation, h eq (4>), where jiy = 0. If the sheets are separated by h < h eq , they will be repelled while if h > h eq , they 
will be attracted. In Fig. 3b we plot h eq as a function of phase difference and we see that the equilibrium separation 
decreases monotonically with decreasing <p and that when the sheets are in phase the vertical force acting on them is 
strictly attractive. Indeed, for <p = 0, in the limit h — > we see that f 2y — (37r/4)(Dei — De2)/(1 + De 2 )[25]. 



Coupled dynamics 

In the idealized two-dimensional case studied here, the swimmer mobility in the vertical direction is strictly zero 
and hence only motion in the horizontal direction occurs. In the slender-body limit, which is the one relevant for the 
dynamics of three-dimensional flagellar filaments of swimming cells, the viscous mobility in the direction perpendicular 
to the length of the flageUum is about half of that in the parallel direction [23]. In order to propose a simple model 
for the coupled vertical/horizontal motion of the sheets, wc proceed to use this ratio in our model, and simply assume 

Given the behavior of the vertical force, we expect the swimmers to be pushed apart slightly if their phase difference 
is large, then, as the phase difference decreases, to be attracted to a final synchronized conformation where the flagella 
are as close together as possible, as seen experimentally [6]. We can numerically integrate both differential equations, 
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FIG. 4: Coupled time evolution of the distance between the swimmers, h(t) (a), and the phase difference, <f>{t) (b). The initial 
condition is h = 2, and the mobility in y is taken to be half of that in x: 4>o = 37r/4 (dotted line), 4>o = t/2 (dashed line), 
(f> = tt/8 (solid line). With e = 0.1, Dei = 100 and De 2 = 10. 

(29) and (37), to obtain the coupled time evolution of h(t), shown in Fig. 4a, and <j>(t), shown in Fig. 4b, for an initial 
separation h = 2. We see that for a small enough initial angle, the sheets are monotonically attracted to each-other 
(0o = 7r/8, solid line). However, for larger initial phase differences, the sheets are initially repelled, before reaching 
a maximum separation, and eventually being drawn together closely. This is illustrated for </>o = tt/2 (dashed line) 
and 4>o = 37r/4 (dotted line). The time scale for the evolution of the phase angle is similar to the constant separation 
case, and all initial conformations converge to the stable in-phase conformation. 

CONCLUSION 

In this paper we used a two-dimensional model to analytically address the synchronization of two swimmers in a 
viscoelastic (Oldroyd-B) fluid. In Newtonian fluids, a front-back asymmetry in the swimmer's waveform is required 
for synchronization. In contrast, in a viscoelastic fluid, phase- locking occurs even for swimmers displaying front /back 
symmetry. The two swimmers are driven into a stable in-phase conformation where a minimum of mechanical energy 
is dissipated, contrary to the Newtonian case where the stable conformation can be either in-phase or opposite-phase 
depending only on the waveform geometry. In addition, the evolution to a phase-locked conformation in a viscoelastic 
fluid occurs on asymptotically faster time scales than in a Newtonian fluid. 

From a biological standpoint, the results of our model indicate that, for example, mammalian spermatozoa pro- 
gressing through cervical mucus would be expected to synchronize passively, thereby reducing the work they are doing 
against the surrounding fluid as compared to when swimming isolated. This net energy savings could then potentially 
be used to increase their wave speeds, and hence swimming speed, as is observed experimentally [6]. The asymp- 
totically larger forces between swimmers might also lead to large-scale coherence in the dynamics of cell suspensions 
which is more pronounced in complex fluids than in Newtonian environments. 
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The other symmetry, with respect to the :r-axis, is always assumed to be true in order to enforce swimming along a straight 
line. 

h is allowed to decrease to zero only when = 0, as otherwise the sheets would overlap. 



